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Abstract 

The  ability  to  reliably  detect  coronary  artery  disease  based  on  the 
acoustic  noises  produced  by  a  stenosis  can  provide  a  simple,  non-invasive 
technique.  Current  research  exploits  the  shear  wave  fields  in  body  tissue 
to  detect  and  analyze  coronary  stenosis.  A  mathematical  model  of  this 
system  couples  the  generation  of  these  acoustic  noises  with  the  propaga¬ 
tion  of  the  sound  and  shear  waves  through  the  chest  cavity.  In  our  initial 
investigations  we  consider  a  one-dimensional,  homogeneous  viscoelastic 
model.  A  quasi-linear  viscoelastic  stress-strain  relationship  was  proposed 
by  Fung  [8]  for  a  variety  of  biological  tissues.  Though  an  effective  model, 
this  formulation  presents  significant  computational  difficulties  in  dynamic 
situations.  We  present  several  alternate  constitutive  relations,  based  on  an 
internal  variable  formulation,  that  approximate  Fung’s  constitutive  rela¬ 
tion  well  when  optimized.  More  importantly,  results  from  the  correspond¬ 
ing  dynamic  models  match  well  with  simulated  data  of  wave  propagation 
through  a  homogeneous  soft  tissue-like  gel. 

Mathematics  Subject  Classification:  65M32,  74D10,  74J30 

Keywords  and  Phrases:  Propagating  shear  waves,  nonlinear  viscoelastic 
materials,  inverse  problems. 

1  Introduction 

Coronary  artery  disease  (CAD)  is  the  buildup  of  plaque  (cholesterol,  calcium, 
and  platelets)  in  the  inner  wall  of  coronary  arteries  known  as  stenosis  in  which 

’Center  for  Research  in  Scientific  Computation,  Box  8205,  North  Carolina  State  University, 
Raleigh,  NC  27695-8205 

tMedAcoustics  Inc, 5540  Centreview  Drive,  Suite  114,  Raleigh,  NC  27606 
tMedAcoustics  Inc, 5540  Centreview  Drive,  Suite  114,  Raleigh,  NC  27606 
§  Center  for  Research  in  Scientific  Computation,  Box  8205,  North  Carolina  State  University, 
Raleigh,  NC  27695-8205 

'Center  for  Research  in  Scientific  Computation,  Box  8205,  North  Carolina  State  University, 
Raleigh,  NC  27695-8205 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

15  AUG  2000 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2000  to  00-00-2000 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Modeling  and  Computation  of  Propagating  Waves  from  Coronary 

_ _ 

5b.  GRANT  NUMBER 

LHCllUSCb 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation,  Raleigh,  NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

21 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Acoustic  Waves 


Turbulence  Stenosis 


Figure  1:  Turbulent  blood  flow  generated  by  a  stenosis. 


blood  flow  is  restricted  and  the  oxygen  supply  to  the  heart  muscle  is  decreased. 
The  end  result  of  arterial  stenosis  is  permanent  damage  to  the  heart  muscle.  It 
is  estimated  that  CAD  affects  more  than  12  million  people  in  the  United  States, 
and  it  is  the  single  leading  cause  of  death  and  premature  permanent  disability 
[2].  As  such,  the  detection  and  treatment  of  CAD  is  a  high  priority. 

It  is  well  known  that  arterial  stenoses  produce  sounds  due  to  turbulent  blood 
flow  in  partially  occluded  arteries  (see  e.g.,  [1]  or  [27]).  In  principle,  turbulent 
normal  wall  forces  exist  at  and  downstream  from  arterial  stenosis.  These  wall 
forces,  which  are  extremely  small,  exert  a  pressure  on  the  wall  of  the  artery, 
which  then  causes  a  small  displacement  in  the  surrounding  soft  tissue  (see  Fig¬ 
ure  1).  The  vibrations  of  the  surrounding  body  tissues,  which  occur  in  two 
forms,  a  compressional  wave  and  a  shear  wave,  produce  sounds  [23],  [27].  In 
larger  arteries  such  as  the  carotid  arteries,  these  acoustic  sounds  can  be  de¬ 
tected  by  physicians  using  a  stethoscope.  However,  detecting  acoustic  signals  in 
smaller  arteries  deep  inside  the  body  has  proved  difficult  for  two  reasons:  these 
acoustic  noises  attenuate  significantly  as  they  travel  through  the  intervening  tis¬ 
sues  [23],  and,  moreover,  many  complex  sounds  within  the  body  can  overwhelm 
conventional  acoustic  detection  systems. 

Current  detection  techniques  include  the  angiogram  which  is  a  reliable,  yet 
expensive,  invasive  technique  and  prone  to  interobserver  variability  (see  e.g.,  [9], 
[17]).  Ultra- fast  CT  techniques  are  also  employed;  this  is  a  non-invasive  imaging 
technique  effective  in  detecting  and  scoring  the  severity  of  calcium  deposits  in 
the  coronary  arteries.  CT  testing  equipment  is  very  expensive.  However,  the 
biggest  limitation  is  that  the  method  only  detects  calcium  deposits  and  not  the 
soft  plaques  that  make  up  many  of  the  most  dangerous  lesions. 

The  ability  to  reliably  detect  coronary  artery  disease  based  on  the  acous¬ 
tic  noises  can  provide  a  simple,  non-invasive  approach  [18].  Current  research 
exploits  the  shear  wave  fields  in  body  tissue  to  detect  and  analyze  coronary 
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stenosis  [15].  At  low  frequencies  (<2  kHz),  shear  wave  attenuation  and  wave 
propagation  speed  is  low,  which  proves  useful  for  studying  propagation  char¬ 
acteristics  and  estimating  mechanical  properties  of  tissue-like  media.  Recent 
techniques  measure  the  shear  wave  propagation  at  the  surface  of  the  chest  wall 
using  a  multiple  array  of  sensors.  Phased  array  signal  processing  is  then  used  to 
image  the  underlying  volume  and  locate  stenosis  sites  within  the  coronary  tree. 

A  complete  mathematical  model  for  this  acoustic  wave  system  couples  two 
separate  processes:  (1)  the  generation  of  the  sound  and  shear  waves  at  the 
arterial  stenosis  and  transmission  through  the  arterial  wall,  and  (2)  the  prop¬ 
agation  of  the  sound,  or  shear  waves,  through  the  chest  cavity  to  the  external 
acoustic  sensors.  Modeling  the  wave  propagation  through  the  chest  cavity  is  a 
nontrivial  task.  It  is  well  accepted  that  body  soft  tissue  medium  behaves  like 
a  viscoelastic  medium  [8].  In  addition,  the  acoustic  waves  will  travel  through 
ribs  and  at  least  two  sorts  of  tissue:  lung  tissue  and  muscular  connective  tissue. 
These  spatial  inhomogeneities  in  the  propagating  medium  cause  spatial  changes 
in  the  medium’s  propagation  speed.  These  variations  induce  refractive  effects, 
which  bend  the  rays  along  which  the  wave  propagates  and  ultimately  leads  to 
multi-paths. 

In  this  paper  we  address  the  second  process:  the  shear  wave  propagation 
through  a  viscoelastic,  heterogeneous  medium.  To  start,  however,  we  first  an¬ 
alyze  the  most  simple  case,  i.e.,  a  one-dimensional,  homogeneous  viscoelastic 
model  for  which  experimental  data  is  readily  available.  The  corresponding 
physical  model  is  depicted  in  Figure  2;  the  synthetic  gel  has  material  prop¬ 
erties  similar  to  those  of  soft  tissues  and  the  tube  mimics  an  arterial  vein  with 
stenosis. 

The  motion  of  a  viscoelastic  body  is  governed  by  the  laws  of  conservation 
of  mass  and  momentum,  the  stress-strain  (constitutive)  relations,  plus  bound¬ 
ary  and  initial  conditions.  A  quasi-linear  viscoelastic  stress-strain  relationship 
which  captures  well  the  viscoelastic  characteristics  for  a  variety  of  biological 
tissues  was  proposed  by  Fung  [8]  for  soft  tissues.  The  constitutive  relation  in¬ 
corporates  a  continuous  spectrum  of  relaxation  times  into  the  model,  and  gives 
the  stress  in  terms  of  a  Boltzmann  integral  which  depends  on  the  history  of  the 
elastic  response  rate.  Though  effective  in  modeling  quasi-static  creep  and  relax¬ 
ation,  this  formulation  presents  significant  computational  difficulties  in  dynamic 
situations. 

In  an  effort  to  avoid  some  of  these  computational  challenges,  we  investigate 
the  effectiveness  of  using  alternative  stress-strain  formulations  to  model  time 
dependent  shear  wave  propagation  through  a  homogeneous  medium.  We  will 
present  several  simple  formulations  that,  when  optimized,  approximate  well 
Fung’s  constitutive  relation.  More  importantly,  results  from  the  corresponding 
dynamic  models  match  well  with  simulated  data  of  wave  propagation  through 
a  homogeneous  soft  tissue-like  gel. 

The  organization  of  the  paper  is  as  follows.  We  first  develop  the  model  equa¬ 
tions  in  Section  2.  In  Section  3  we  introduce  alternative  stress-strain  constitu- 
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Figure  2:  The  ID  homogeneous  viscoelastic  model. 


tive  relations  and  compare  with  Fung’s  formulation.  The  physical  and  dynamic 
models  are  given  in  Section  4,  and  we  discuss  the  computational  method  for 
integrating  the  models  in  Section  5.  In  Section  6  we  present  model  simulations 
and  compare  with  data,  and  follow  with  conclusions  in  Section  7. 


2  Model  Formulation 


We  first  develop  the  one  dimensional  model  equations  based  on  Newton’s  law  of 
motion,  the  principle  of  conservation  of  mass,  and  laws  of  thermodynamics.  In 
the  following,  let  u.  a,  and  X  denote  the  Cartesian  shear  displacement,  shear 
stress,  and  body  force  per  unit  volume,  respectively.  Let  p  denote  the  mass 
density  and  v  denote  the  shear  velocity  du/dt.  X  is  a  source  term  for  the 
acoustic  sound  propagation. 

From  the  conservation  of  mass,  we  have  the  continuity  equation 

I  +  =  <*> 

Conservation  of  momentum  implies  that 


dpv 

dt 


d  2\  v 

+  s*'"  > =  s  +  x 


Using  equation  (1)  in  the  left  side  of  the  above  equation,  we  have  the  Eulerian 


4 


equation  of  motion  of  a  continuum 


dv  dv 

m  +vfa 


da 

dx 


+  X. 


Since  we  are  interested  in  only  small  displacements,  we  ignore  the  higher  order 
term  pvdv/dx  to  obtain 

d2u  da 

pW  =  d^  +  x-  (2) 

We  next  obtain  a  constitutive  equation  that  relates  stress  and  strain,  and 
eventually  stress  and  displacement.  The  medium  of  interest  is  soft  tissue  (arter¬ 
ies,  muscle,  skin,  lung,  etc...)  which  is  considered  viscoelastic  material.  That  is, 
soft  tissue  exhibits  the  characteristics  of  creep  (a  slow  progressive  deformation  of 
a  material  under  constant  stress),  relaxation  (a  gradual  decrease  of  stress  when 
the  material  is  held  at  constant  deformation),  and  hysteresis  (the  stress-strain 
curve  in  the  loading  process  is  usually  different  than  in  the  unloading  process). 

The  quasi-linear  relation  proposed  by  Fung  [8]  describes  the  quasi-static 
behavior  of  many  biological  soft  tissues  (see,  e.g.,  [5], [10], [11],  [19], [20], [21], [25], 
[28], [29], [30]).  This  relation  between  stress  and  strain  in  simple  elongation  is 
given  by 


,  /■  d<re(A(r))  , 

alt)  =  Git-  t)  ,  dr. 
Jo  dr 


(3) 


where  G(t)  and  ae( A(-))  are  the  reduced  relaxation  function  and  the  elastic 
response,  respectively.  Here,  the  relaxation  function  is  assumed  to  be  separable 
in  time  and  strain. 

By  definition,  the  elastic  response  is  the  tensile  stress  instantaneously  gen¬ 
erated  in  the  tissue  when  a  step  function  of  stretching  A  is  imposed  on  the 
specimen.  For  many  materials,  it  may  be  approximated  by 


^-  =  a(ae+fi),  a  e  ( 1 )  —  0 , 


(4) 


where  a  and  /3  are  constants  to  be  estimated  from  data.  For  a  justification  of 
this  approximation,  see  [8],  p.  279.  We  recall  [3],  [8]  that  the  stretch  ratio  A  is 
related  to  the  infinitesimal  strain  e  (and  displacement  u)  by  the  equation 

A2  =  1  +  2^=  (l  +  ^)  =(l+e)2. 


Solving  equation  (4)  for  ae,  we  then  have 

ae(X(t))  =  -p  +  pea9u/9x. 

The  reduced  relaxation  function  proposed  by  Fung  is  given  by 


G(t) 


-{ 


1  +  C 


El(±)-El(-) 
?2  Ti 


n 

l+cln 

J  J 

Vn/J 

1  -1 


(5) 


(6) 
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where  Ei  (z)  is  the  exponential  integral  function  defined  by 


Ei 


dt. 


The  parameters  C,  n ,  and  T2  are  constants  determined  from  data;  the  constant 
C  represents  the  degree  to  which  viscous  effects  are  present,  n  and  r2  represent 
the  fast  and  slow  viscous  time  phenomena,  respectively  [22].  Fung’s  formulation 
provides  a  continuous  spectrum  of  relaxation.  This  is  a  departure  from  other 
models  that  commonly  express  the  relaxation  function  as  the  sum  of  exponential 
functions  where  each  exponent  is  identified  with  the  rate  constant  of  a  relaxation 
mechanism,  i.e., 


G(t) 


E  Cte-** 
EC 


According  to  Fung,  such  a  model  will  have  discrete  relaxation  rate  constants 
and  corresponds  to  a  discrete  hysteresis  spectrum;  this  is  incompatible  with  the 
observation  that,  for  many  biological  soft  tissues,  the  hysteresis  loop  is  nearly 
independent  of  the  strain  rate  within  several  decades  of  the  rate  variation. 

Equations  (2),  (3),  (5),  and  (6)  comprise  an  integro-differential  system.  In 
summary,  the  one-dimensional  equations  for  wave  propagation  through  a  ho¬ 
mogeneous,  viscoelastic  medium  are  the  equation  of  motion  (2)  and  Fung’s 
stress-strain  relation  (3)  with  the  elastic  response  and  reduced  relaxation  func¬ 
tion  defined  by  (5)  and  (6),  respectively.  Because  the  constitutive  equation  (3) 
is  in  integral  form  with  a  complicated  kernel  described  by  equation  (6) ,  it  must 
be  approximated  computationally.  This  requires  the  storage  of  the  history  of 
the  kernel  function  G(t)  in  small  time  increments.  Furthermore,  this  time  dis¬ 
cretization  might  be  considerably  smaller  than  the  time  discretization  used  in 
the  approximation  of  the  equation  of  motion  (2)  (so  that  the  errors  in  computing 
the  constitutive  equation  do  not  contribute  to  the  errors  in  approximating  the 
acoustic  energy).  This  is  especially  true  in  the  case  where  the  elastic  response 
<re  is  highly  oscillatory,  as  in  our  study  here.  For  computational  efficiency,  we 
propose  alternatives  to  Fung’s  constitutive  equation  in  the  following  section. 


3  Constitutive  Law  Formulations  and  Compar¬ 
isons 

In  this  section  we  propose  three  computationally  tractable  alternatives  to  Fung’s 
quasi-linear  constitutive  equation  for  soft  tissues  and  compare  them  numerically 
to  Fung’s  formulation.  The  idea  is  to  consider  internal  variable  models  (Maxwell 
solids  in  differential  form)  as  in  [12]  and,  more  recently,  in  investigations  of 
nonlinearities  and  hysteresis  arising  from  a  class  of  magnetorheological-based 
smart  elastomers  [3],  [4].  The  basic  assumption  is  that  one  has  a  finite  number 
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of  “internal  strain”  variables,  =  1  driven  by  the  infinitesimal 

strain  e(t)  =  X (t)  —  1.  In  this  formulation  the  stress  is  given  by 

N 

<7(*)  =  5Ze^)'  (7) 

i—1 

The  internal  strains  can  be  modeled  dynamically  as 

^  =  -¥i+Ci^(A(f))1  j  =  1, . . . ,  N.  (8) 

Using  the  variation  of  constants  formula,  solutions  to  the  above  system  of  dif¬ 
ferential  equations  with  zero  initial  conditions,  ej( 0)  =  0,  are  given  by 

e>(t)  =  j'cje-'A-Q^aeMsVds,  (9) 

Note  that  the  internal  variable  approach  (7)  is  equivalent  to  the  integral  for¬ 
mulation  (3)  with  an  exponential  form  for  the  reduced  relaxation  function,  i.e., 
G(t)  =  J2^=i  Cje~Vit.  The  difference  in  these  two  formulations  is  only  in  their 
numerical  implementation:  using  (3)  one  deals  with  numerical  approximation 
of  an  integral  function  without  any  consideration  of  internal  dynamics,  whereas 
using  the  internal  variable  approach  (7)  one  has  to  numerically  solve  a  decou¬ 
pled  system  of  differential  equations.  Computationally,  one  expects  that  the 
internal  variable  model  approach  is  much  more  efficient  in  terms  of  both  speed 
and  data  storage. 

More  generally,  the  internal  strains  might  be  modeled  by  nonlinear  dynamics 
of  the  form 

^  + c^ixm  j  =  v  ■  ,n.  (io) 

However,  this  formulation  cannot  be  written  in  Boltzmann  form. 

We  will  compare  Fung’s  constitutive  relation  (3)  with  several  constitutive 
relations  resulting  from  the  internal  variable  approach.  Let  N  =  1  in  (7),  then 
o{t)  =  ei(f),  G{t)  m  Cie~vit,  and  a  is  obtained  as  the  solution  of  the  following 
linear  ordinary  differential  equation  (ODE) 

^  =  -Vla(t)  a(  0)=0.  (11) 

The  constants  C±  and  V\  are  unknown  constants  to  be  determined  later;  ae( X) 
is  given  in  (5).  Similarly,  if  we  let  N  =  2,  then  a(t)  =  (t)  +  ti (t)  and  G(t)  ~ 

Cie~Vlt  +  The  internal  strain  variables  ei,  ti  satisfy  the  following 

linear  ODEs,  respectively, 

+  Ci-^-(X(t)),  ei(0)=0  (12) 

^  =  -V2t2{t)+C^{  X(t)),  e2(0)=0. 
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The  parameters  C\,  C2,  Vi,  and  v 2  are  the  corresponding  unknown  constants 
to  be  determined  later.  Considering  the  nonlinear  dynamics  formulation  (10) 
with  N  =  1,  we  have  in  a  third  case  a(t)  =  (t)  and  a  is  the  solution  of  the 

following  nonlinear  ODE 

^=g(a(t))+K^(  A(t)),  cr(0)  =  0.  (13) 

We  assume  that  g(z )  =  1  ai^i(z)  is  a  piecewise  linear  function  of  a.  In  this 

case  there  are  five  constants  to  be  determined:  K.  a±,  02,  03,  04. 

To  determine  the  effectiveness  of  the  internal  variable  formulations,  we  com¬ 
pared  the  stress  obtained  numerically  from  each  relation  given  above:  Fung’s 
stress-strain  integral  equation  (3)  using  kernel  (6),  the  ODE  from  using  a  one- 
exponential  kernel  (11),  from  using  a  two-exponential  kernel  (12),  and  the  piece- 
wise  linear  ODE  (13).  In  each  case,  we  set  ux(x,t)  =  —  msin(u)t),  with  to  =  0.1 
and  u)  £  [100, 600],  the  frequency  range  of  interest.  Values  for  a  and  /3  were  set 
as  a  =  0.1  and  f3  =  1.0  (see  [6]).  To  compute  Fung’s  kernel  we  set  c  =  0.05, 
7~i  =  0.005,  and  72  =  50  (again,  see  [6],  [22]).  We  will  assume,  for  the  sake  of 
comparison,  that  the  constitutive  equation  from  Fung  is  exact. 

Each  of  the  equations  was  solved  numerically  in  MATLAB.  For  the  inte¬ 
gral  expression  we  used  an  adaptive  recursive  Newton  Cotes  8  panel  rule.  The 
ODEs  were  solved  using  the  solver  ODE15s,  a  variable  order  method  for  solv¬ 
ing  stiff  ODEs.  We  initially  set  the  unknown  constants  to  values  such  that  the 
corresponding  kernels  were  decent  approximations  of  Fung’s  kernel  (6).  The 
constants  were  then  optimized  using  the  Nelder-Mead  simplex  algorithm  [14]  so 
as  to  best  approximate  the  a  obtained  using  Fung’s  kernel. 

Figure  3  depicts  the  stresses  at(t)  obtained  by  each  of  the  four  constitutive 
equations  with  frequency  cu  =  600Hz.  Since  the  function  is  very  oscillatory, 
we  show  Ui(t )  only  for  t  £  [0.09,0.1]s.  Figure  4  is  the  same  plot,  but  over 
a  smaller  time  interval  to  indicate  more  clearly  the  numerical  differences.  In 
each  plot,  a  obtained  from  Fung’s  method  is  indicated  with  a  solid  line  labeled 
Vfung,  the  one-exponential  method  is  indicated  with  a  dashed  line  (aexpi),  the 
two-exponential  method  is  indicated  with  a  dashed-dot  line  (aexp 2),  and  the 
nonlinear  method  is  indicated  with  a  dotted  line  (ani). 

The  differences  between  a  using  Fung’s  kernel  and  the  alternate  constitutive 
equations  with  optimized  constants  are  slight.  We  emphasize  that  the  unknown 
constants  were  optimized  to  match  a,  and  not  so  that  the  exponential  kernels 
best  match  Fung’s  kernel.  In  general,  more  exponential  terms  in  the  kernel  are 
needed  to  obtain  good  agreement  with  Fung’s  kernel. 

We  conclude  from  this  numerical  experiment  that  the  alternate  constitutive 
equations  match  quite  well  with  Fung’s  stress-strain  relation.  Since  the  alternate 
relations  are  in  differential  form,  the  formulation  is  easy  to  implement  and  no 
longer  requires  the  cumbersome  storage  of  kernel  history  data.  In  the  following 
sections  we  include  the  various  stress-strain  equations  in  the  dynamic  model 
and  compare  results  with  simulated  data.  First,  we  present  the  physical  and 


8 


Figure  5:  The  simplified  physical  model. 


dynamic  models,  together  with  boundary  and  initial  conditions,  and  discuss  the 
numerical  method  used  to  integrate  the  models. 

4  Physical  and  Dynamic  Models 

The  simplified  physical  model  is  a  cylindrical  gel  mold  (such  as  those  used  in 
physical  experiments  at  MedAcoustics)  with  a  surgical  tube  passing  through 
the  center  axisymmetrically.  A  source  disturbance  is  generated  within  the  tube. 
Shear  waves  propagate  through  the  gel,  and  the  shear  displacement  is  measured 
(by  sensors)  at  the  gel  outer  surface.  The  gel  possesses  the  material  character¬ 
istics  of  soft  tissue.  The  disturbance  within  the  tube  represents  a  source  due  to 
stenosis.  Let  R\  and  R2  be  the  inner  and  outer  radius  of  the  gel,  respectively, 
(see  Figure  5).  We  assume  that  the  source  disturbance  is  time  dependent,  has 
only  a  radial  component  (i.e.,  it  is  axisymmetric) ,  and  the  disturbance  is  a  pure 
shear  force.  The  outer  surface  of  the  gel  is  a  free  surface.  The  gel  is  initially  at 
rest. 

The  one  dimensional  equations  governing  the  shear  stress  (a (t,  r))  and  shear 
displacement  (u(t,  r ))  are 


d2u 

da 

PW 

dr 

a(t) 

=  fcit 

Jo 

~  s)^-it3ea9ru)  ds 

OS 

(14) 

(t,Ri) 

=  fit), 

a{t,R2 )  =  0, 
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where  G(t)  is  Fung’s  reduced  relaxation  function,  j3eadvU  represents  the  elastic 
response,  and  f(t)  is  the  input  shear  stress.  Since  we  are  interested  in  the 
feasibility  of  a  more  computationally  efficient  stress-strain  relation,  we  will  not 
integrate  the  dynamic  model  with  Fung’s  stress-strain  relation. 

Instead,  in  light  of  the  simulation  results  for  alternative  stress  formulations  in 
the  previous  section,  we  concentrate  on  the  internal  strain  variable  constitutive 
relations.  Using  equation  (11)  as  the  stress-strain  relation  (i.e.,  G(t)  «  Cie-"1*), 
we  find  that  system  (14)  is  equivalent  to 

d2u  _  da 

PW  = 

f  =  -^  +  ^4^.)  (15) 

a(t,Ri)  =  fit),  a(t,R2)=  0, 


The  unknown  constants  are  C±,  v\,  and  a;  since  C±  and  (3  always  appear  as  a 
product,  we  set  9  =  1  and  let  C\  be  the  unknown  constant.  Using  equation 
(12)  as  the  stress-strain  relation  ( G(t )  «  Cie~Vlt  +  C2e~V2t),  system  (14)  is 
equivalent  to 


d2u  d  . 

<’T»  ~  HR1  +  ei) 

^  =  -W2  +  C, (16) 

ei(t,l?i)  =  7  fit),  eiit,R2)  =  0, 

^it,Ri)  =  (1-7  )f(t),  e2(t,R2)=0, 


where  ait,r)  =  ei(t,r)  +  e2 (t,r),  and  7  is  a  free  parameter  used  to  set  the 
boundary  conditions.  The  unknown  constants  in  this  system  are  C\,  C2,  v\, 
v2 ,  and  a.  Finally,  if  we  use  the  piecewise  linear  constitutive  relation,  equation 
(13),  then  we  have  the  following  nonlinear  system 


d2u  _  da 

PW  ~ 

=  »w+Ai<e“8-“)  <17) 

ait,Ri)  =  fit),  a{t,R2)=  0, 

where  gia)  =  as  before.  The  unknown  constants  are  K,  ax,  a2, 

az,  04,  and  a. 

In  the  next  section  we  describe  the  numerical  method  used  to  integrate  each 
of  the  approximate  systems. 
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5  Computational  Methods 


We  use  the  same  method  to  numerically  integrate  each  of  the  approximate 
systems  presented  above.  For  simplicity,  we  present  the  method  with  respect 
to  the  one-exponential  system  given  in  (15);  the  two-exponential  and  piecewise 
linear  systems  are  treated  similarly. 

First,  we  rewrite  the  ( u,a )  system  as  a  first  order  system  in  v,  w,  and  p 
where  v  =  du/dt ,  w  =  du/dr,  and  p  =  a  -  C(5eaw  (or  pi  =  e,  —  Ci/3eaw).  The 
corresponding  (v,w,p)  system  for  system  (15)  is 

-viitt  +  CxPe™) 

^O.  +  Ci/Je-)  (18) 

dv 
dr’ 

which  is  now  in  conservative  form  and  is  essentially  hyperbolic. 

Physical  boundary  conditions  are  prescribed  in  terms  of  a,  but  are  translated 
into  boundary  conditions  for  w  by  analytically  integrating  the  stress-strain  equa¬ 
tion  in  system  (15)  over  (0 ,t),  and  solving  for  w(t,r)  at  the  boundaries  (r  =  f?i 
and  r  =  R2).  The  boundary  conditions  for  w  are  then 


dp 

dt 

dv 

dt 

dw 

dt 


w{t,Ri) 

w(t,R2) 


1 


In 


'fit)  -  /( 0)  +  +  Vl  f*  f(s)  ds' 


a  y 

w(0  ,R2). 


C i/3 


Since  the  gel  is  initially  at  rest,  the  initial  conditions  are  v(0,  r)  =  0,  u>(0,  r)  =  0, 
and  /i(0 ,r)  =  —Ci/3. 

To  numerically  integrate  the  system,  subject  to  the  above  boundary  and 
initial  conditions,  we  employ  a  MacCormack  finite  difference  scheme  which  is 
second-order  in  space  and  time.  Boundary  conditions  for  v  and  //  are  updated 
using  local  direction  cosines  at  the  domain  boundaries  (see,  e.g.,  [13]  or  [24] 
on  hyperbolic  systems).  If  the  local  direction  points  into  the  domain,  then  the 
update  uses  the  new  boundary  data  (i.e.,  w(tn+1,Ri)  or  w(tn+1,  R2))-  Oth¬ 
erwise,  the  local  direction  points  out  of  the  domain,  and  boundary  values  are 
updated  using  data  from  the  domain  interior  and  second-order  differences  to 
approximate  any  spatial  derivatives. 

The  stability  of  the  method  is  verified  through  a  Von  Neumann  analysis 
of  the  linearized  ( fx,v,w )  system  about  the  initial  condition.  We  apply  the 
MacCormack  scheme  to  the  linearized  model  and  then  substitute  in  planar  waves 
of  the  form  V  =  yel(kr+wt) .  w  =  w^(kr+ut)  _  ancj  ^  kr+wt )  _  Le^.  x  =  eiwAt 

be  the  amplification  factor,  and  let  0  =  j Ar,  r  =  At/Ar  and  7  =  Cf3a/p. 
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Nontrivial  solutions  for  p ,v,w  imply  the  following  condition  for  x> 


(x-i) 


(x  -  l)2  +  (x 


-«{ 


vAt(l  —  — y-)  +  2r27(l  —  cos  9) 


+T'7  {"a,<1  “  ^)(1 " cos8)  +  (t27^)(1  “ cos9)2 


=  0. 


Solving  this  expression  for  x,  we  have  x  =  1,  or 


X 


1  —  t27(1  cos#) 


vA  t 


(1 


vA  t. 


v2A t2  vA t  v2A t2  i/Af  .  2 

±  \l — - — (1 - Y - - — r27(l  -  cos#r  -  t27(1 - rsin  9 

4  2  4  2 


The  scheme  is  stable  when  |x|  <  1.  Since  x  depends  on  seven  parameters: 
C,  v,  a,  j3,  p,  At/Ar,  and  At,  we  easily  check  the  stability  of  the  scheme  for  a 
particular  parameter  set  using  the  above  expression  for  x,  rather  than  determine 
a  domain  of  stability.  For  the  computational  experiments  described  in  the  next 
section,  we  verified  that  the  parameter  set  produces  a  stable  scheme. 

We  performed  a  grid  refinement  to  verify  the  second-order  accuracy  of  the 
scheme.  This  was  done  with  the  forced  (p,v,w)  system  by  choosing  an  exact 
solution  and  setting  the  forcing  functions  appropriately.  We  chose  the  following 
exact  solutions  for  u  and  a, 


u(t,r)  =  —m(  1  —  cos  (ut)) 


Ri 


i?2  —  Rl 


r(t,  r)  =  msin(cjt) 


Tl> 


—  Rl 


where  m  is  the  oscillation  amplitude,  and  to  is  the  oscillation  frequency.  For 
these  calculations,  m,  =  0.1  and  uj  €  [100,600]  Hz.  The  exact  solutions  for  the 
{p.  v,  w )  system  are  calculated  as  p  =  a  —  C/3eaUr ,  v  =  ut,  and  w  =  ur. 

To  perform  the  grid  refinement,  we  fix  the  ratio  ££  and  repeatedly  halve 
the  time  step  and  the  spatial  increment.  The  truncation  error  for  each  variable 
is  computed  as  the  L norm  of  the  difference  between  the  exact  and  numerical 
solutions.  For  example,  the  truncation  error  for  p  is 


En{p)  =  ||  p(T,n)  —  pf  ||oo 


where  /z(T,  r* )  and  pf  are  the  exact  and  numerical  solutions,  respectively,  at 
the  point  (T,  rY).  The  superscript  N  =  T/  At  indicates  the  number  of  numerical 
time  increments,  and  r*  =  R\  +  iAr,  i  =  1, . . . ,  Nr.  represents  the  spatial  grid 
points.  The  subscript  n  indicates  the  level  of  refinement;  initially  n  =  0.  With 
each  refinement,  we  increment  n  and  double  the  total  number  of  time  steps 
(N)  and  the  number  of  grid  points  (Nr).  Setting  ^=l/20sm-1.,  T  =  5s, 
R2  —  Ri  =  lm.,  and  Nr  =  40  for  the  initial  convergence  test,  we  compute  the 
ratio  En(-)/En+ 1(-)  for  each  state  variable.  The  results  for  v  are  given  below. 
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Table  1:  Grid  Refinement  Analysis  for  v 


n 

l|£n(«)||oo 

PnMIlWPn+lMIloo 

IT 

0.1295 

1 

0.0321 

4.03 

2 

0.0080 

4.01 

IT 

0.0020 

4.00 

4 

5.0054e-04 

3.99 

A  ratio  of  4  indicates  the  code  is  second  order  accurate.  Results  for  w,  and 
H  are  similar. 

Hence  the  scheme  is  indeed  stable  and  second-order  accurate.  In  the  next 
section  we  will  compare  this  system,  as  well  as  the  two  exponential  and  piecewise 
linear  systems,  with  published  data. 


6  Model  Simulations 

To  measure  the  ability  of  the  approximate  systems  to  match  experimental  data, 
we  compare  the  results  of  the  dynamic  models  presented  in  Section  4  with 
simulation  data  derived  from  the  experimental  results  of  Verburg.  This  is  a  first 
step  in  the  ultimate  goal  of  matching  shear  displacement  data  measured  on  the 
the  outer  surface  of  the  gel. 

In  [26]  Verburg  presents  wave  speeds  per  frequency  for  shear  waves  propa¬ 
gating  through  a  homogeneous  medium  with  mechanical  properties  similar  to 
the  lung  tissue  located  between  the  heart  valves  and  the  chest  wall.  Output 
data  is  generated  from  the  Verburg  data  in  the  following  way.  First,  the  input 
signal  is  a  sum  of  weighted  sine  waves  with  a  sampling  rate  of  12288  Hz  and 
A /  =  6  Hz.  Using  the  Verburg  wave  speeds,  a  phase  shift  is  calculated  for  each 
input  frequency.  The  output  signal  is  then  the  weighted  sum  of  sine  waves  with 
velocity  phase  delays.  The  resulting  simulation  data  does  not  account  for  atten¬ 
uation  due  to  damping  or  dispersion  loss.  Also,  the  input  and  output  data  sets 
are  the  same  order  derivative  in  time.  Since  the  dynamic  model  expects  input 
shear  stress,  we  interpret  the  simulation  data  as  force  in  (scaled  by  a  constant) 
and  acceleration  out  (see  Figure  6  for  a  graph  of  the  input  function). 

To  compare  the  model  with  this  simulated  data,  we  solve  the  associated  best 
fit  inverse  problem.  In  other  words,  we  make  an  initial  guess  for  the  unknown 
constants  in  the  model,  then  optimize  the  constants  using  the  Nelder-Mead 
algorithm.  The  Nelder-Mead  cost  function  is  the  L2  norm  of  the  difference 
of  the  simulated  output  data  and  the  model  acceleration  data  at  the  outer 
boundary  (r  =  R2).  The  model  acceleration  data  is  generated  by  solving  the 
forward  problem,  i.e,  by  numerically  integrating  the  dynamic  model.  Note,  the 
acceleration  data  is  obtained  as  a  finite  difference  approximation  of  dv/dt. 

In  Figures  7-9  we  compare  the  simulated  Verburg  data  and  the  optimized 
dynamic  model  data.  Each  figure  contains  two  plots.  The  top  plot  compares 
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Figure  6:  The  input  shear  stress  function  a(t,Ri )  =  f(t). 


the  normalized  acceleration  data  (normalized  since  the  simulated  data  does  not 
account  for  dissipation) .  The  bottom  plot  compares  the  normalized  fast  Fourier 
transform  (FFT)  of  the  output  data.  This  plot  indicates  which  frequencies  are 
excited  by  each  model.  In  all  plots,  the  simulated  Verburg  data  is  indicated 
with  a  solid  line;  the  model  data  is  indicated  with  a  dashed  line. 

Figure  7  exhibits  data  from  the  one-exponential  system  (15),  Figure  8  presents 
data  from  the  two-exponential  system  (16),  and  data  from  the  piecewise  linear 
system  (17)  is  shown  in  Figure  9. 

Note  that  the  unknown  constants  are  optimized  only  with  respect  to  the  ac¬ 
celeration  data.  That  is,  the  fit-to-data  uses  the  acceleration  data;  the  FFT  of 
the  optimized  model  is  then  compared  to  the  FFT  of  the  data.  Any  additional 
matching  of  the  FFT  data  is  an  independent  indicator  that  the  alternate  mod¬ 
els  are  reasonable.  In  each  case  the  optimized  systems  achieve  the  correct  time 
lag  in  the  acceleration  data.  None  of  the  model  simulations  contain  the  small 
oscillations  that  appear  for  t  >  0.03s,  probably  because  the  internal  strain  vari¬ 
able  models  contain  dissipation  whereas  the  simulated  data  does  not.  Clearly, 
the  optimized  two- exponential  model  produces  the  best  approximation  with  the 
simulated  data;  differences  in  this  case  are  nearly  indiscernible. 

7  Conclusions 

In  the  case  of  simple  one-dimensional  shear  wave  propagation  through  a  homo¬ 
geneous  soft  tissue  (viscoelastic)  medium,  we  have  demonstrated  that  dynamic 
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simulation  using  an  internal  variable  based  constitutive  relation  matches  well 
with  data.  Moreover,  the  internal  strain  variable  formulations  lead  to  systems  of 
first  order  differential  equations,  rather  than  an  integro-differential  model,  mak¬ 
ing  the  alternate  models  more  computationally  tractable.  It  should  be  noted 
that  these  alternate  constitutive  relations  may  not  match  well  with  Fung’s  re¬ 
lation  (3)  in  all  respects  of  viscoelastic  behavior.  However,  the  main  goal  in 
the  efforts  reported  here  is  to  match  shear  displacement  data  after  propaga¬ 
tion  through  a  soft  tissue- like  medium.  In  this  respect,  the  optimized  models 
with  internal  variable  strain  formulations  are,  we  believe,  adequate.  The  one 
and  two-exponential  models  correspond  to  standard  kernel  formulations  in  vis¬ 
coelastic  theory  (i.e.,  a  finite  series  of  Kelvin  models).  It  is  quite  possible  that 
the  simulated  Verburg  data  does  not  capture  all  nonlinear  behavior  that  might 
be  present  in  complex  experimental  data.  Hence  our  models  with  linear  internal 
dynamics  and  nonlinear  coupling  with  the  infinitesimal  strain  produce  a  good  fit 
to  this  data.  However  we  also  achieve  decent  agreement  between  data  and  the 
model  with  nonlinear  internal  dynamics  (as  represented  by  the  piecewise  linear 
model) ,  which  offers  promise  for  use  of  this  class  of  models  with  more  complex 
data.  Future  efforts  will  include  modeling  of  the  two-dimensional  propagation 
problem.  In  this  case  there  is  a  variety  of  experimental  data  available  to  be 
employed  in  model  validation. 

8  Acknowledgments 

This  research  was  supported  in  part  (H.T.B.  and  S.W.)  by  the  U.S.  Air  Force  of 
Scientific  Research  under  grant  F-49620-98-1-180.  The  authors  wish  to  thank 
Dr.  Z.  Li  and  K.  Coffey  for  many  valuable  discussions. 


References 

[1]  M.  Akay,  Noninvasive  detection  of  coronary  artery  disease  using  advanced 
signal  processing  methods,  PhD.  dissertation,  Rutgers  University,  Piscat- 
away,  NJ  1990. 

[2]  American  Heart  Association,  1999  heart  and  stroke  statistical  update, 
American  Heart  Association,  1998. 

[3]  H.T.  Banks,  G.A.  Pinter,  L.K.  Potter,  M.J.  Gaitens,  and  L.C.  Yanyo, 
’’Modeling  of  quasi-static  and  dynamic  load  responses  of  filled  viscoelas¬ 
tic  materials”,  CRSC-TR98-48,  Raleigh,  NC,  December  1998;  Modeling: 
Case  Studies  from  Industry,  (E.  Cumberbatch  and  A.  Fitt,  eds.),  Cam¬ 
bridge  Univ.  Press.,  to  appear. 

[4]  H.T.  Banks,  G.A.  Pinter,  L.K.  Potter,  B.C.  Munoz,  and  L.C.  Yanyo,  ’’Esti¬ 
mation  and  control  related  issues  in  smart  material  structures  and  fluids” , 


19 


CRSC-TR98-2,  Raleigh,  NC,  January  1998;  Opt.  Techniques  and  Applica¬ 
tions ,  (L.  Caccetta,  et.  al.,  eds.),  Curtain  Univ.  Press,  p.  19-34,  1998. 

[5]  Y.L.  Chen  and  Y.C.  Fung,  Stress-strain  history  relations  of  rabbit  mesen¬ 
tery  in  simple  elongation,  Biomechanical  Symposium,  ADM-2,  ASME , 
1973,  9-10. 

[6]  L.J.M.G.  Dortmans,  A.A.H.J.  Sauren,  and  E.P.M.  Rousseau,  Parameter  es¬ 
timation  using  the  quasi-linear  viscoelastic  model  proposed  by  Fung,  Trans, 
of  the  ASME,  J.  Biomech.  Engr.,  106,  1984,  198-202. 

[7]  Y.C.  Fung,  A  First  Course  in  Continuum  Mechanics  Prentice-Hall,  Inc., 
New  Jersey,  (1977). 

[8]  Y.C.  Fung,  Biomechanics:  Mechanical  Properties  of  Living  Tissues 
Springer- Verlag,  New  York,  (1993). 

[9]  J.E.  Galbraith,  M.L.  Murphy,  and  N.  de  Soyza,  Coronary  angiogram  inter¬ 
pretation.  Interobserver  variability,  Jama,  240,  1978,  2053-2056. 

[10]  R.C.  Haut  and  R.W.  Little,  A  constitutive  equation  for  collagen  fibers,  J. 
Biomechanics,  5,  1972,  423-430. 

[11]  R.  B.  Jenkins  and  R.W.  Little,  A  constitutive  equation  for  parallel-fibered 
elastic  tissue,  J.  Biomechanics,  7,  1974,  397-402. 

[12]  A.R.  Johnson,  A.  Tessler,  and  M.  Dambach,  ’’Dynamics  of  thick  elastic 
beams”,  J.  Engr.  Materials  and  Tech.,  119,  1997,  273-278. 

[13]  J.  Kevorkian,  Parital  Differential  Equations:  Analytical  Solution  Tech¬ 
niques,  Wadsworth,  Inc.,  California,  (1990). 

[14]  J.C.  Lagarias,  J.A.  Reeds,  M.H  Wright,  and  P.E.  Wright,  Convergence 
properties  of  the  Nelder-Mead  simplex  method  in  low  dimensions,  SIAM 
J.  Optim.,  9,  1999,  112-147. 

[15]  MedAcoustics  Inc.,  internal  documentation.  Raleigh,  NC. 

[16]  I.  Nigul  and  U.  Nigul,  On  algorithms  of  evaluation  of  Fung’s  relaxation 
function  parameters,  J.  Biomechanics,  20,  1987,  343-352. 

[17]  S.E.  Nissen  and  J.C.  Gurley,  Application  of  intravascular  ultrsound  for  de¬ 
tection  and  quantitation  of  coronary  atherosclerosis,  Int.  J.  Cardiac  Imag¬ 
ing,  6,  1991,  165-177. 

[18]  N.L.  Owsley  and  A.J.  Hull,  Beamformed  nearfield  imaging  of  a  simulated 
coronary  artery  containing  a  stenosis,  IEEE  Trans.  Med.  Imaging,  17,  1998, 
900-909. 


20 


[19]  J.G.  Pinto  and  P.J.  Patitucci,  Visco-elasticity  of  passive  cardiac  muscle, 
ASME  J.  Biomechanical  Eng.,  102.  1980,  57-61. 

[20]  E.P.M.  Rousseau,  A.A.H.J  Sauren,  M.C.  van  Hout,  and  A. A.  van 
Steenhoven,  Elastic  and  viscoelastic  material  behavior  of  fresh  and 
gluteraldehyde-treated  porcine  aortic  valve  tissue,  J.  Biomechanics,  16, 
1983,  339-348. 

[21]  A.A.H.J.  Sauren,  M.C.  van  Hout,  A. A.  van  Steenhoven,  F.E.  Veldpaus,  and 
J.D.  Janssen,  The  mechanical  properties  of  porcine  aortic  valve  tissues,  J. 
Biomechanics ,  16,  1983,  327-338. 

[22]  A.A.H.J.  Sauren  and  E.P.M.  Rousseau,  A  concise  sensitivity  analysis  of  the 
quasi-linear  viscoelastic  model  proposed  by  Fung,  ASME  J.  Biomechanical 
Eng.,  105,  (1983),  92-95. 

[23]  J.L.  Semmlow,  W.  Welkowitz,  J.  Kostis  and  J.W.  MacKenzie,  Coronary 
artery  disease  correlates  between  diastolic  auditory  characteristics  and 
coronary  artery  stenoses,  IEEE  Trans.  Biomed.  Eng.,  BME-30,  1983, 136- 
139. 

[24]  J.C.  Strickwerda,  Finite  Difference  Schemes  and  Partial  Differential  Equa¬ 
tions,  Wadsworth,  Inc.,  California,  (1989). 

[25]  T.T.  Tanaka  and  Y.C.  Fung,  Elastic  and  inelastic  properties  of  the  canine 
aorta  and  their  variation  along  the  aortic  tree,  J.  Biomechanics,  7,  1974, 
357-370. 

[26]  J.  Verburg,  Transmission  of  vibrations  of  the  heart  to  the  chestwall,  Adv. 
Cardiovasc.  Phys.  5,  1983,  84-103. 

[27]  J.  Verburg  and  E.  van  Vollenhoven,  Phonocardiography:  Physical  and  tech¬ 
nical  aspect  and  clinical  uses,  in  Noninvasive  Physiological  Measurements 
(ed.  P.  Rolfe),  Academic  Press,  London,  (1979). 

[28]  S.L.Y.  Woo,  B.R.  Simon,  S.C.  Huei,  and  W.H.  Akeson,  Quasi-linear  visco¬ 
elastic  properties  of  normal  articular  cartilage,  ASME  J.  Biomechanical 
Eng.,  102,  1980,  85-90. 

[29]  S.L.Y.  Woo,  M.A.  Gomez,  and  W.H.  Akeson,  The  time  and  history- 
dependent  viscoelastic  properties  of  the  canine  medial  collateral  ligament, 
ASME  J.  Biomechanical  Eng.,  103,  1981,  293-298. 

[30]  S.L.Y.  Woo,  Mechanical  properties  of  tendons  and  ligaments:  I.  Quasi¬ 
static  and  nonlinear  viscoelastic  properties,  Biorheology,  19,  1982,  385-396. 


21 


